Two-Component Nonlinear Schrodinger Models with a Double- Well Potential 
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We introduce a model motivated by studies of Bose-Einstein condensates (BECs) trapped in 
double-well potentials. We assume that a mixture of two hyperfine states of the same atomic 
species is loaded in such a trap. The analysis is focused on symmetry-breaking bifurcations in the 
system, starting at the linear limit and gradually increasing the nonlinearity. Depending on values 
of the chemical potentials of the two species, we find numerous states, as well as symmetry-breaking 
bifurcations, in addition to those known in the single-component setting. These branches, which 
include all relevant stationary solutions of the problem, are predicted analytically by means of a two- 
QQ \ mode approximation, and confirmed numerically. For unstable branches, outcomes of the instability 

, development are explored in direct simulations. 

o : 

^ . I. INTRODUCTION 

' The nonlinear Schrodinger (NLS) equation is a ubiquitous partial differential equation (PDE) with a broad spectrum 
of applications, including nonlinear optics in the temporal and spatial domains, matter waves, plasmas, water waves, 
' and some biophysical models [l|, 0| • One of the most fundamental applications of the NLS equation stems from its 
direct relevance as an accurate mean-field model (known as the Gross-Pitaevskii equation, GPE, in that context) 
of Bose-Einstein Condensates (BECs) 0, The GPE usually includes an external potential accounting for the 
magnetic, optical or combined confinement of dilute alkali vapors that constitute the BEG [H]- Basic types of such 
trapping potentials include parabolic and spatially periodic ones (the latter is created, as an optical lattice, by the 
^ interference of counter-propagating laser beams). The NLS equations with similar potentials are also relevant models 
^ for optical beams in graded- index waveguides and periodic waveguiding arrays ^gj, lS] • 

^ , A setting that has recently drawn much interest in the context of BEGs is based on a double- well potential (DWP), 

\ which originates from a combination of the two above-mentioned types of potentials. It was realized experimentally in 
^ . [Si] , leading to particularly interesting phenomena including tunneling and Josephson oscillations (for a small number 
^ ' of atoms), or macroscopic quantum self-trapping, with an asymmetric division of atoms between the two wells, for 
, a large number of atoms. Numerous theoretical studies of DWP settings have been performed in parallel to the 
• experimental work [1, [l^, [Hi [H, [H, [3, HI, [H, [13 • They addressed finite- mode reductions, analytical results for 
specially designed shapes of the potential, quantum effects, and other aspects of the theory (in particular, tunne ling 
between vortex and antivortex states in BEG trapped in a two-dimensional (2D) anisotropic harmonic potential [l8| 
Jlf^ \ belongs to the latter category). 

An interesting generalization of the DWP is provided by multi-well potentials. In particular, nonequilibrium BEG 
states and generation of quantum entanglement in such settings were recently studied in detail theoretically in . 

Also considered were 2D [2^ [2l| and 3D extensions of DWP settings, which add one or two extra dimensions to 
the model, either without any additional potential, or with an independent optical lattice acting in these directions. 
The extended geometry makes it possible to consider solitons, self-trapped in the extra dimension(s), which therefore 
5J] ■ emerge as effectively ID j^^, [2l| or 2D [l^ dual-core states. The solitons may be symmetric and antisymmetric, 
\ as well as ones breaking their symmetry between the wells through bifurcations. These states may be described by 
simplified systems of linearly coupled ID |2fli] or 2D [20, 22] PDEs, or, in a more accurate form, effectively ID solitons 
can be found as solutions to the full two-dimensional PDE, that takes into regard a particular form of the DWP 
(which depends on the transverse coordinate and is uniform in the longitudinal direction, allowing the solitons to 
self-trap in that free direction) [2l|. 

DWPs are also relevant to nonlinear-optics settings, such as the twin-core self-guided laser beams in Kerr media 
[23} , optically induced dual-core waveguiding structures in a photorefractive crystal ^] , and DWP configurations for 
trapped light beams in a structured annular core of an optical fiber [2^. As concerns DWPs with extra dimensions, 
their counterparts in fiber optics are twin-core fibers [26| and fiber Bragg gratings [27'| that were shown to support 



symmetric and asymmetric solitons. In addition, also investigated was the symmetry breaking of solitons in models of 
twin-core optical waveguides with non-Kerr nonlinear terms, viz, quadratic |28j and cubic-quintic (GQ) p9| . All these 
optical model are based on systems of linearly coupled ID PDEs. Also in the context of nonlinear optics, a relevant 
model is based on a system of linearly coupled complex Ginzburg-Landau equations with the GQ nonlinearity, that 
gives rise to stable dissipative solitons with broken symmetry [s^l • 

In addition to the linearly coupled systems of two nonlinear PDEs motivated by the DWP transverse structures, 
several models of triangular configurations, which admit their own specific modes of the symmetry breaking, have also 
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been introduced in optics, each model based on a system of three nonhnear PDEs with symmetric hnear couphngs 
between them. These include tri-core nonlinear fibers [3l| and fiber Bragg gratings l32i] . as well as a system of three 
coupled Ginzburg-Landau equations with the intra-core nonlinearity of the CQ type [33| . 

The significant interest in the DWP settings has also motivated the appearance of rigorous mathematical results re- 
garding the symmetry-breaking bifurcations and the stability of nonlinear stationary states ^34] . A rigorous treatment 
was also developed for the dynamical evolution in such settings (ssl . [sH] • 

Our objective in this work is to extend the realm of DWPs to multi-component settings. This is of particular 
relevance to experimental realizations of BEG - e.g.,, in a mixture of different hyperfine states of ®^Rb [3, 4] (see 
also very recent experiments reported in Ref. [37| and references therein). This extension is relevant to optics too, 
where multi-component dynamics can be, for instance, realized in photorefractive crystals [ssj . In the one-component 
setting, the analysis of the weakly nonlinear regime has given considerable insight into the emergence of asymmetric 
branches from symmetric or anti-symmetric ones, in the models with self-focusing and self-defocusing nonlinearity, 
respectively), and the destabilization of the "parent" branches through the respective symmetry-breaking bifurcations, 
as well as the dynamics initiated by the destabilization 0, [13, 113, Ea, [3 113] ■ 

In the present work, we aim to extend the analysis of the DWP setting to two-component systems. In addition to 
the finite-mode approximation, which can be justified rigorously ji4] under conditions that remain valid in the present 
case, we use numerical methods to follow the parametric evolution of solution branches that emerge from bifurcations in 
the two-component system. We observe that, in addition to the branches inherited from the single-component model, 
new branches appear with the growth of the nonlinearity. We dub these new solutions "combined" ones, as they 
involve both components. In fact, the new branches connect some of the single-component ones. Furthermore, these 
new branches may undergo their own symmetry-breaking bifurcations. The solutions depend on chemical potentials 
of the two species, fii and ^2 (in terms of BEG), and, accordingly, loci of various bifurcations will be identified in 
plane (Ati,/Lt2). 

The paper is structured as follows. In Section II, we present the framework of the two-component problem, 
including the two-mode reduction (in each component), that allows us to considerably simplify the existence and 
stability problems. The formulation includes also a physically possible linear coupling between the components, but 
the actual analysis is performed without the linear coupling. In section III, we report numerical results concerning 
the existence and stability of the new states, as well as the evolution of unstable ones. In section IV, we summarize 
the findings and discuss directions for further studies. 



II. ANALYTICAL CONSIDERATION 

As a prototypical model that is relevant both to BEGs [1, [3| and optics [6| , we take the normalized two-component 
NLS equation of the following form: 

iut — Lu + Kv + (T(|up -l- (7|wp)u — Hiu, 
ivt = Lv + KU + cr(|wp + g\u\'^)v — 112V, 

where u and v are the wave functions of the two BEG components, or local amplitudes of the two optical modes, 
/ii.2 are chemical potentials in BEG or propagation constants in the optical setting, k is the coefficient of the linear 
coupling between the components, and 

L ^-{1/2)81 + V{x) (2) 

is the usual single-particle operator with trapping potential V{x). The repulsive or attractive character of the non- 
linearity in BEG (self-defocusing or self- focusing, in terms of optics) is defined by cr = -t-1 and cr = — 1, respectively, 
while ag is the coefficient of the inter-species interactions in BEG, or cross-phase modulation (XPM) in optics. 

In the case of (for instance) two hyperfine states with \F,niF >= > and |2,1 > in ^^Rb, both coefficients 

of the self-interaction and the cross-interaction coefficient are very close, although their slight difference is critical 



in accounting for the immiscibility of the two species (see, e.g.,, 37[ and references therein). However, for problems 
considered below, the difference does not play a significant role, hence we set 5 = 1 in Eqs. ([1]), which corresponds to 
the Manakov's system, i.e., one which is integrable, in the absence of the potential, V{x) = in Eq. ^ [aH^. In 
fact, V{x) is the DWP, which we compose of a parabolic trap (with corresponding frequency fl) and a sech^-shaped 
barrier (of strength Vq and width w) in the middle: 

V{x) = (l/2)f7V + Fosech^ (x/W) . (3) 
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The calculations presented below have been performed for Q = OA, Vq = I and W = 0.5, in which case the first 
two eigenvalues of linear operator ^ are ujq = 0.1328 and uji — 0.1557. However, we have checked that the generic 
picture presented below is insensitive to specific details of the potential, provided that it is symmetric. 

The linear coupling in Eqs. ([T]), via the terms proportional to k, accounts for the possibility of intcrconvcrsion 
between the two hyperfine states in the BEC mixture, which may be induced, e.g.,, by appropriate two-photon pulses 
in the situation considered in Ref. [37l |. In optics, the linear coupling is induced by a twist of the optical waveguide, 
if the two wave components correspond to orthogonal linear polarizations of light. 

To define an appropriate basis for the finite-mode expansion of solutions to Eqs. ([1]), we first consider the following 
eigenvalue problem. 



(4) 



We denote the ground and first excited eigenstates of the matrix operator in Eq. which appertain to the two 
lowest eigenvalues, as {uo,vo} and {uiVi} (these two states are, as usual, spatially even and odd ones, respectively). 
The Hermitian nature of the linear operator in Eq. Q ensures that the eigenvalues are real, and the eigenfunctions 
can also be chosen in a real form. In fact, in the case that we will examine here, we will have uq = vq and ui = vi, 
but the different notations will be kept for the u and v components, to demonstrate how the nonlinear analysis can 
be carried out in the general case. Then, the two-mode approach is based on the assumption that, in the weakly 
nonlinear case, the solution is decomposed as a linear combination of these eigenfunctions, i.e., 

u{x,t) = co{t)uo{x) + ci{t)ui{x), 

V{x,t) = C2{t)vo{x) + C3{t)vi{x). 

This assumption can be rigorously substantiated, with an estimate for the accuracy, upon imposing suitable conditions 
on the smallness of the nonlinearity, and properties of the rest of the linear spectrum [34]. Therefore, as shown in 
[3^ . such an approach can be controUably accurate for small N — (|itp + jvp) dx, although below we examine 
its comparison with numerical results even for values of of 0(1). 

It is relevant to mention that an alternative approach to the decomposition of the two-component wave functions 
could be based on using the basis of symmetric and antisymmetric wave functions, uq ± vq and ui ± wi. However, 
the resort to this basis does not make the final equations really simpler. On the other hand, unlike the approach 
developed below, the use of the alternative basis would make it much harder to estimate the error of the finite-mode 
approximation, and thus rigorously substantiate the approximation, cf . Ref. [s^l . 

Substituting ansatz ([5]) in Eqs. ((T]), and projecting the resulting equations onto eigenfunctions {uq, vq} and {uiVi}, 
one can derive at the following nonlinear ODEs for the temporal evolution of the complex amplitudes of the decom- 
position, co,i,2,3: 



icQ = {uq - /ii)co + KC2 + o-poooolcol^co + rooii(ciCQ + 2co|ci|^)] 

-I- erg [Eoooo Co I C2p + rooii(co|c3p -I- C1C2C3 + cic^cj)], 

ici = {ui - fii)ci + KC3 + o-piiiilcipci + rooii(cocJ + 2|copci)] 

+ o-.g[riiiiCi|c3|^ + rooii(ci|c2|^ + C0C2C3 + C0C2C3)], 

iC2 = (wo - M2)C2 + KCO + O-[r0000|c2pC2 + rooil(c3C2 + 2c2|c3p)] 

-I- crg[roooo|copC2 4- rooii(|cipC2 -I- cocjca -I- cjjcica)], 

zc's = (t^i - ^2)C3 + KCi + 0-[riiii|c3pC3 + rooil(c2C3 + 2|c2pC3)] 

+ o-.g[riiii|cipC3 -I- Tooiidcopcs + coclc2 + C0C1C2)]. 



(6) 



(7) 



(8) 



(9) 



In these equations, Tijki = Ui{x)uj{x)uk{x)ui{x)dx are the so-called matrix elements of the nonlinear four-wave 
interactions, with Tijki = when i+j+k+l is odd. 

Seeking for real fixed points to Eqs. (O - ([5]), aj{t) = pj, we reduce the equations to an algebraic system, 

MiPo = (t^oPo + KP2) + ^(rooooyOo + SrooiiPoP?) ^^^^ 
+ crg(ToooQPQP2 + Tooiipopj + 2rooiipip2P3), 
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MiPi = i^iPi + i^Pa) + '^(riiiiPi + SrooiiPoPi) ^^^^ 
+ ag{Tiiiipipl + TqqiipipI + 2rooiipoP2P3), 

P2P2 = {l^0P2 + Kpo) + CT(rooooP2 + 3rooilP2P3) ^^^^ 
+ o-ff(roooo/0o/02 + ^ooiipIp2 + 2rooilPoPl/03), 

f^2P3 = (wiP3 + Kpi) + aiTiiiipl + SToohpIps) 

/ 2 2 \ \^) 

+ crg[Tiiiip^p3 + TqquPqPs + 2rooiiPoPiP2), 

The simplest solution to the above equations can be obtained in the form of p2 = po and p3 = pi, for pi = p2- With 
this substitution, the remaining equations are 

pi- uja = a-(rooooPo + Srooiipf ), ^^^^ 
Pi-LOi ^ d-{Tiiiipl + 3rooiiPo)> 

where Cjq = ujq + k, loi = uji + k, a = a{l + g). System ([14]) can be solved as a linear one for pQ i. For given values 
of r and Lu, a physical range of chemical potential pi is that which provides physical solutions for Pq i. 

In what follows, we focus on the case of k = 0, since k ^ imposes the condition pi = P2, while we are interested 
in more general solutions. With k = 0, chemical potentials pi and p2 may be different, allowing us to explore the 
two-parameter solution space {pi, P2), for both signs of the nonlinearity, self-attractive (cr = —1) and the self-repulsive 
('^ = +!,)• 

The simplest possible branches of solutions are single-mode ones, with only one of amplitudes pj different from 
zero. Branches of these solutions can be easily found from Eqs. (fTO)) - (IT3)) . 
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The first one of these branches contains the projection only onto the even (symmetricO eigenfunction of the first 
component, and will hereafter be accordingly denoted SI. Similarly, the second branch has a projection onto the odd 
(antisymmetric) mode of the first component, and will be named ANl. The third and fourth solutions are similar 
modes in the second component, to be denoted S2 and AN2, respectively. 

In addition to these modes, there exist asymmetric states in each of the components (similar to the one-component 
model [li,|33). 



2 i- iiii\hi-i - ^0 - 'J'- ooii\h>-i - ^1 n 
Po = ?7F2 — ^ ' P2 = 0. (16) 



„2 _ riii i(^i - lvq) - 3rooii(/ii - uji) 
o'(riiiiroooo — 9r§oii) 

2 _ TooooImi ^ ^1) ^ 3rooii(/^i — Wo) 



Pi = -r^ f Xf^^ > P3 = 0. (17) 

iiiii 0000 — oOlli 

Since these solutions do not exist in the linear limit, they have to bifurcate at a non-zero value of the amplitude, from 
either symmetric or antisymmetric branches (jlSp . In fact, such solutions bifurcate, depending on the sign of a and 
coefficients T, either from SI at 

(cr) 3Tooii{iLJo — UJi) 

Pi =^1--^ ^ , (18) 

i 0000 — 'J'- 0011 



or from ANl at 



(cr) 3rooil(wi — Wo) 

Pi -^0--^ ^ ■ (19) 

i 1111 — oi 0011 
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The asymmetric solutions given by Eqs. ((T6)) and (fT7|) will be denoted ASl. Similarly, there is an asymmetric branch 
of solutions in the second component, hereafter denoted AS2, with 



2 rini(^2 ^o) - 3roon(/i2 - uji) 
P2 - -7f f ^f^^ ' Po = (20) 

iiiii 0000 ^ I 



0011-' 



2 rooOo(A^2 — iJJl) — 3rooil(^2 — ^o) n /oi \ 

Pz = -7F ?^7^2 ^ ' (21) 



c(riiiiroooo ^ 



0011^ 



and conclusions that can be deduced about its bifurcation are similar to those concerning ASl. 

In addition to these obvious solutions, there exist others, that can be obtained as solutions to Eqs. pO| - (fT3l) and 
involve both components (obviously, they have no counterparts in the single-component model). They will be denoted 
as "combined" ones (CI, C2, C3 and C4) in what follows. Although they cannot be expressed by simple analytical 
formulas, unlike the S, AN and AS branches that were defined above, they can be easily found as numerical solutions 
of Eqs. (fT0)l -p3 l) . and will be compared to full numerical results in the following section. 



III. NUMERICAL RESULTS 



We start the numerical analysis by examining the self-focusing case. Because of the complexity of the bifurcation 
diagrams in the plane of the chemical potentials of the two components, (/xi,/i2), we will illustrate the bifurcation 
phenomenology by showing cross-sections of the full two-parameter bifurcation diagram for varying /ii, at different 
fixed values of /i2. 

In connection to the single-component states presented above, four distinct scenarios have been found, when the 
second chemical potential, /i2, appears. The first (most complex) scenario takes place at /i2 < ^^2^^ (with some 
critical value ^^^2'^^), when branch AS2 exists (in addition to S2 and AN2). The second case is observed in region 
/Xg^"^^ < ^2 < ^Oi in which case only S2 and AN2 exist (recall Wq and Wi are two lowest eigenvalues of operator ^). 
The third scenario is found in the interval of cjq < ^^2 < (^i, where only branch AN2 may exist, and, finally, the fourth 
one takes place at loi < fi2, where there are no single-component solutions for the second field (in the case of the 
self- attractive nonlinearity) . 

We describe the most complex of these scenarios in Fig. [T]and subsequently summarize differences from this case 
in the remaining three simpler regimes. In Fig. [1] the top panel presents the full numerically generated bifurcation 
diagram, in terms of total norm A'^ = -I- jwp) dx (which is proportional to the total number of atoms in 

both species in the BEC model, or the total beam power in optics) as a function of ni for fixed /i2 = 0.1. The 
horizontal lines in this diagram represent single-component solutions AS2, S2 and AN2 in the second field (in order 
of the increasing norm), which are, obviously, unaffected by the variation of /zi of the other component, that remains 
"empty" in these solutions. As the asymmetric branch AS2 exists at this value of fj.2 (in fact, at all values in region 
ft2 < Mo'^'^'* ) 1 it is stable (as the one generated by the symmetry-breaking bifurcation in the single-component equation 
[is Hi)) while branch S2 is destabilized by the same bifurcation. In addition to these branches represented solely by 
the second component, with the variation of /ii we observe the emergence of the single-component branches, SI and 
ANl, in the first field, from the corresponding linear limit, at ni = loq and /xi = coi. In addition, past critical point 
fi\ = 0.1214, the symmetry-breaking-induced branch ASl arises in the first component, destabilizing its symmetric 
"parent branch" SI, in agreement with the known results for the one-component model. 

An additional explanation is necessary about the special case of fii = 112 (= 0.1 in the case of Fig. [T|). In that case, 
we observe that branches SI and S2, as well as ANl and AN2, and also ASl and AS2, collide. At this value of the 
parameter, stationary equations ([1]) (with ut = Vt = 0) are degenerate (recall we consider the case of g — 0), making 
any solution in the form of au = bv is possible, including pairs (a, b) — (0, 1) (which corresponds to the branches 
considered above that contain only the first component), (a, 6) — (1,0) (i.e., the second-component branches), and 
(a, b) = (1, 1) (where the two components are identical). Therefore, in this degenerate case, the points of collisions of 
the branches in the bifurcation diagram of Fig. [factually denote a one-parameter family of solutions, rather than an 
isolated solution. 

The key new feature of this bifurcation diagram in comparison to its one-component counterpart is the presence 
of the four new branches of solutions, CI through C4. Two of them are found to connect ("bridge") two single- 
component branches of the solutions belonging to different components, and two other species of the new solutions 
bifurcate from the bridging ones. More specifically, as can be seen in the bottom left panel of Fig. [T] (which is a 
blowup of the one in the top left corner), branch C2 interpolates between antisymmetric branch AN2 of the second 
component and symmetric branch SI of the first component (hence, one may naturally expect that the solutions 
belonging to C2 have symmetric and antisymmetric profiles in the first and second components, respectively, see Figs. 
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[5][7] below. Furthermore, similarly to the bifurcation of the asymmetric branches from the symmetric ones in the 
single-component models, there is such a bifurcation which occurs to the C2 branch, destabilizing it and giving rise 
to a stable branch CI of combined solutions, with asymmetric profiles instead of the symmetric and antisymmetric 
waveforms in the first and second components of C2, respectively (see Fig. [S] below). 

Naturally, there is a similar pair of two-component branches emerging due to the bridging of S2 and ANl. As 
expected, the two-component branch C4, which is responsible for the bridging per se, features symmetric and anti- 
symmetric profiles in its second and first components. Finally, it is clearly seen in the bottom right panel of Fig. [1] 
(which is a blowup of another relevant fragment of the top left panel) that a symmetry-breaking bifurcation occurs 
on branch C4, which destabilizes it and gives rise to branch C3 of stable asymmetric two-component solutions. It is 
worthy to note that this emerging asymmetric branch C3 terminates through a collision with the single-component 
asymmetric branch AS2, as can also be seen in the top left panel of Fig. [TJ 

Similar to the emergence of ASl from SI, the bifurcation of CI from C2 and of C3 from C4 are of the pitchfork type. 
Therefore, we actually have two asymmetric branches of types CI and C3, which are mirror images to each other 
(accordingly, the additional asymmetric waveforms can be obtained from those displayed in Fig. [7] by the reflection 
around x — 0). 

Comparing this complex bifurcation picture with predictions of algebraic equations (|10p - (|13p . we conclude that the 
picture is captured in its entirety by the two-mode approximation, as shown in the top right panel of Fig. [T] Although 
some details may differ (notice, in particular, slight disparity in the scales of the two diagrams in the top panels of Fig. 
[T]), the overall structure of the diagrams is fully reproduced by the algebraic equations. As a characteristic indicating 
the accuracy of the approximation, it is relevant to compare the bifurcation points. In the finite-mode picture, the 
bifurcation of ASl from SI occurs at fii = 0.1211, CI bifurcates from C2 at = 0.0780, and the bifurcation of C3 
from C4 happens at /ii = 0.1225. The full numerical results yield the following values for the same bifurcation points: 
/ii = 0.1214, fii = 0.0781 and fii = 0.1224, respectively, with the relative error < 0.25%. 

Two additional descriptions of these bifurcations, which help to understand the emergence of the combined 
branches, are presented in Figs. [5] and [3] The former one shows norms = \u\^dx and iV^, — \v\'^dx, 

and the latter figure shows the total number of atoms in the left and right wells, Nl — jl^ + I^^P) dx and 
= /q°° (I^P -I- dx, as functions of the chemical potential ^i. The left panels are produced by the numerical 
solution of the stationary version of underlying GPEs, Eqs. H]), while the right panels correspond to the finite-mode 
approximation. Notice a sharp (near vertical) form of the combined branches C2 and C4 in the former figure, which 
indicates a short interval of their existence, and their change of stability upon the collision with (or, more appropri- 
ately, the bifurcation of) CI and C3, respectively. In the latter figure, it is worthwhile to notice the similarity of the 
curves for Nl and Nn appertaining to branches SI, ANl, S2 and AN2, in contrast with the differences between the 
respective curves for ASl, AS2 and combined branches CI and C3, which indicates the asymmetry between the states 
trapped in the two wells in the latter case. 

(cr) 

Having addressed the most complex bifurcation scenario observed at /i2 < ^2 > now turn to the three remaining 

cases, namely, < ^J'2 < ^o, i^o < 1^2 < '^i, and ivi < /i2- Full numerically generated bifurcation diagrams for 

each of these cases can be found in the left panels of Fig. 01 while the corresponding results produced by algebraic 
equations (|10p - (|13p are presented in the right panels. Once again, we notice a very good agreement between the two 

sets of the results. In interval /Xj'^''' < fJ-2 < i-^o, the main difference from the case shown in Fig. [T]is that, as concerns 
the single-component branches belonging to the second field, the asymmetric branch AS2 has not bifurcated from the 
symmetric one S2, therefore S2 is a stable branch. As a result, the two-component branches C4 and C3 do not exist 
in this case (in particular, the existence of C3 is not possible topologically, since it would destabilize C4 which would 
subsequently have to merge with stable branch S2). 

Nevertheless, the behavior of branches CI and C2 remains the same as before. In the middle panels of Fig. 21 
which represents the case of < /i2 < , the situation is simpler in that branch S2 does not exist (in this regime, 
only branch AN2 exists in the second component), hence the presence of C3 and C4, that would connect S2 to ANl, 
is impossible. Finally, in the case of uji < fJ-2, there are no single-component branches in the second field; as a result, 
even the two-component branch C2, joining AN2 and SI, has to disappear (hence CI bifurcating from it cannot exist 
either), and we are therefore left solely with the single-component solutions for u (only the branches SI, ANl and 
ASl are present). 

Examples of solutions representing all the branches considered above are shown in Figs. OITl Figure O represents 
the three branches that have only the u component, namely SI, ANl and ASl in the left, middle and right panels of 
the figure, for /ii = 0.04 and /i2 = 0.1. Notice that, since solution ASl exists for this parameter sets, solution SI is 
unstable, while ANl and ASl are stable. 

The results of the linear-stability analysis around the solutions are shown in the bottom panels of the figure. For 
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FIG. 1: (Color online) Top panel: The norm (corresponding to the number of atoms in BEC and total power in optics) of the 
numerically found solutions of Eq. ([1} (left) and their counterparts predicted by the two-mode approximation (right) for the 
self-attractive nonlinearity {a — —1), as a function of fii, for jj,2 = 0.1. Here and in other figures, the (blue) solid lines and 
(red) dashed ones depict stable and unstable solutions, respectively. The bottom panels are blowups of segments in the top 
left panel where bifurcations involving the new combined solutions occur. SI, ANl, ASl and S2, AN2, AS2 mark, respectively, 
symmetric, antisymmetric, and asymmetric single-component solutions belonging to fields u or v. Symbols CI, C2, C3, C4 
mark branches of the new combined (two-component) solutions, which are defined in the text. 



the purpose of this analysis, perturbed versions of stationary solutions {uq{x),vo{x)} are taken as 

u{x,t) = uoix) + e (c/i(a;)e^* + C/2 (a;)e^**) , 
v{x,t) = vo{x) + e(Vi{x)e^' + V^{x)e^*'^ , 



(22) 
(23) 



where e is an infinitesimal amplitude of perturbations, and the resulting linearized equations for eigenvalue A and 
eigenvector {Ui, U2, Vi, ¥2)'^ are solved numerically. As usual, instability is manifested by the existence of eigenvalue(s) 
A with a non-zero real part. The stability results are shown in terms of the spectral plane (A^, Ai) for A = A,- + iXi, 
hence the instability corresponds to the presence of an eigenvalue in the right-hand half-plane (for example, for branch 
SI in Fig. [5]). Very similar results to those displayed in Fig. [5] can be obtained for solutions S2, AN2 and AS2, in 
which only the second component is nonzero (because they are direct counterparts to those in Fig. [51 they are not 
shown here). 

The combined solutions emerging due to the bridging of the single-component branches are shown in Figs. [6] and 
[71 The former figure shows two instances of branch C2 (right and middle panels), before and after the bifurcation of 
branch CI (the latter one is shown in the left panel). The stability features of these solutions are displayed in the 
bottom panels of the figure. Similar features are shown in Fig. [3 for two-components branches C4 and C3. 

We now turn to the self-defocusing case (corresponding to repulsive interatomic interactions in BEC), with (7 = 1 
in Eqs. (|T]). The results are shown in Figs. [51[51 In this case, we do not show profiles of solutions belonging to various 
branches, as they are very similar to those obtained in the model with the self-attraction, the only difference being 
that their spatial size is larger, due to the self-repulsive nature of the nonlinearity. 
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FIG. 2; (Color online) The norm of the numerically found (left) and approximate two-mode (right) wave functions u (top) and 
V (bottom) of solutions to the stationary version of Eq. ([1} with the attractive nonlinearity (cr = — 1), as a function of /ii, for 
/i2 = 0.10. The notation is the same as in Fig. [1] 



( cr ) 

We again distinguish the main regimes, namely, < Ai2, when all three branches S2, AN2 and AS2 exist; 

wi < /i2 < Ij!'2'^\ for which only S2 and AN2 exist; loq < ^2 < ^i, when only S2 is present; and /Lt2 < Wq, for which 
there is no branch of single-component solutions in the v field. The first (and most complex) of these regimes is 
shown in Fig. [51 in direct analogy to Fig. [T] for the focusing case. Once again, we observe the presence (in addition 
to the three single-component branches in the u field, and three such branches in the v field) of four families of the 
combined solutions. In Fig. [51 C2 (with symmetric u and anti-symmetric v components) bridges stable branch S2 
and unstable one ANl, while branch CI of stable two-component asymmetric solutions bifurcates from C2, making 
it unstable. There is also branch C4 which links SI and AN2, and C3, which bifurcates from C4 and terminates 
by merging into AS2. Details of this picture are shown in the bottom panels of Fig[8l To highlight the accuracy 
of the approximation based on algebraic equations ([TU)) - ([T^ . we display the bifurcation diagram predicted by this 
approximation in the top right panel of Fig. [51 which confirms that all the bifurcations and branches are captured 
within the two-mode framework. As a measure of the agreement between the approximate and numerical results, we 
again present coordinates of the bifurcation points: ASl emerges from ANl at /ii = 0.1684 according to the numerical 
results, while the approximation predicts this to happen at /.ti = 0.1682; further, the bifurcations of CI from C2 and 
C3 from C4 are found numerically to occur at //i = 0.2045 and jii = 0.1561, respectively, while Eqs. pO)l - (fT3|) 
predict the corresponding values fii = 0.2045 and /ii = 0.1561. 

The remaining three regimes and the bifurcation scenarios predicted by the finite-mode approximation are shown 
in the left and right panels of Fig. [9] The top panel corresponds to the case of ui < ^2 < l^2^\ when the absence 
of asymmetric branch AS2 (upon merging into which, C3 would terminate) prevents the existence of the pair of 
two-component branches, C3 and C4, while CI and C2 persist in this case. In the middle panel, corresponding to 
< M2 < ^^r, we observe that only S2 survives among the single-mode branches in field v, and none of the combined 
solutions is present. Finally, as might be naturally expected, no branches with u — Q exist at ^2 < ^0, hence no 
bifurcations of two-component branches are possible (i.e., we get back to the one-component picture). In all these 
three cases, we observe good agreement of the numerically found scenarios with the bifurcation diagrams predicted 
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FIG. 3: (Color online) The norms of the right (top) and left (bottom) parts of the wave functions, Nr = f^°° |itp + \v\'^dx, 
and Nl = f'^^ + \v\'^dx, as obtained from the numerical (left) and approximate two-mode (right) solutions of Eq. ([1} with 
the attractive nonlinearity (cr = — 1), as a function of /ii, for ~ 0.10. The notation is the same as in Fig. [T] Notice that, as 
in the top left panel of Fig. [T] the branches C2 and C4 are actually indiscernible. 



by the two-mode approximation, which are displayed in the right panels. 

An interesting difference between the self- focusing and defocusing nonlinearities is that, as seen in both Figs. [8] and 
[9l in the latter case the bifurcating branches (such as ASl and CI) may become unstable due to a Hamiltonian Hopf 
bifurcation and the resulting oscillatory instability associated with a quartet of complex eigenvalues. This instability, 
which was known in the single-component setting (see, e.g.,, Ref. fl£l|), affects a short dashed (red-colored) interval 
within branches ASl and CI in Figs. [8]and[9l An example of such an unstable configuration and the associated 
spectral plane of the stability eigenvalues are shown in Fig. [Tni 

Finally, we examine the dynamical evolution of the unstable two-component solutions that were revealed by the 
above analysis. The evolution of a typical unstable solution belonging to branches C2 (left panels) and C4 (right 
panels) past the bifurcation points (at which branches CI and C3 emerge, respectively) is shown for the cases of 
the self-attractive and self-repulsivc nonlinearity in Figs. [TT] and 1121 The main dynamical feature apparent in the 
instability evolution is the growth of the asymmetry of the wave functions between the two potential wells. In each 
case, this leads to different wells trapping more atoms (or more power, in terms of optics) in each one of the two 
components. For instance, in the left panel of Fig. [TTl we observe, at i « 300, that the first component features a 
larger norm in the right well, while the second component - in the left one. Due to the Hamiltonian nature of the 
system, the two components do not settle down into a static asymmetric configuration, but rather oscillate between 
the two mirror-image asymmetric states, which (as stationary solutions) are generated by the pitchfork bifurcation. 
In particular, the instability of branch C2 causes oscillations around the two asymmetric states belonging to branches 
CI, and, similarly, the evolution of C4 leads to oscillations around the two states of type C3. 
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FIG. 4: (Color online) The norm of the numerical (left) and approximate two-mode (right) stationary solutions of Eq. ^ for 
the attractive nonlinearity {a = —1) as a function of ni for fi2 = 0.13 (top), /i2 = 0.14 (middle), and ~ 0.16 (bottom). The 
notation is the same as in Fig. [T] 



IV. CONCLUSIONS 



In this work, we have presented the phenomenology and full bifurcation analysis of two-component mixtures trapped 
in the DWP (double- well potential). The model is of straightforward interest to BEC and nonlinear optics (in the 
spatial domain). In our analytical considerations we have developed a two- mode (in terms of each component) 
approximation, that reduces the search for stationary states to solving a set of algebraic equations. The bifurcation 
diagrams obtained in this approximation, which involve all relevant solutions, have been verified by the comparison 
with their counterparts produced by the numerical solution of the full PDE model. The novel feature of the two- 
component setting, in comparison with its previously explored single-component counterpart, is the existence of 
numerous branches of the "combined" solutions that involve both components. These branches emerge from and 
merge into previously known single-component ones. The new branches may combine a symmetric field profile in 
one component and an anti-symmetric one in the other. In addition, asymmetric (in both components) combined 
branches have been found too; they emerge from the symmetric/anti-symmetric two-component states via pitchfork 
bifurcations, similar to how bifurcations of the same type give rise to asymmetric solutions in single-component models. 
The stability analysis of all the considered branches confirms expectations suggested by the general bifurcation theory, 
according to which the pitchfork destabilizes the previously stable "parent" branch, from which the two new ones 
(mirror images of each other) emerge. Direct numerical simulations of unstable symmetric two-component solutions 
(past the bifurcation point) indicate that the instability leads (quite naturally) to oscillations around the asymmetric 
profiles emerging from the bifurcation. 

This study can be extended in several directions. On the one hand, it would be interesting to address this two- 
component setting using the phase-space analysis, similarly to how it was done in [§] . Another possibility is to develop 
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FIG. 5: (Color online) Top panel: profiles of wave function u corresponding to the symmetric (left), antisymmetric (middle) 
and asymmetric (right) branches of the single-component solutions, SI, ANl, and ASl, respectively, in Fig. [T]for /ii = 0.04. 
Bottom panel: the result of the linear-stability analysis around SI (left), ANl (middle) and ASl (right), in the complex plane 
of (ReA, ImA). The existence of an eigenvalue with a positive real part implies instability of the solution. 



an extension of the one-component shooting method of [40| to find an entire set of stationary solutions of the system; 
this is similar to what has been done for two-component optical lattices in the interesting, very recent work of [4l| (also, 
many of the relevant solutions appear quite similar to the ones obtained herein). It should nevertheless be mentioned, in 
connection to 41], that in the present problem there is a straightforward linear limit whose eigenfunctions/eigenvalues 
allow us to construct all the possible solutions that emerge in the presence of nonlinearity. It may also be interesting 
to consider the two-component double-well system in the presence of a spatially modulated nonlinearity in the spirit 
of Ref. [i^l; note that both our semi-analytical and numerical approach could be directly adapted to that setting. 
Finally, it would be particularly interesting to extend this analysis to a four-well, two-dimensional setting (which may 
naturally emerge from a combination of a magnetic trap with a square 2D optical lattice in a "pancake-shaped" BEC). 
Numerous additional solutions and a much richer phenomenology may be expected in the latter case. Some of these 
directions are currently under study. 
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FIG. 7: Profiles of wave functions u and v (top panel's solid and dash-dotted line, respectively) and stability eigenvalues 
(bottom) corresponding to two-component solutions C3 and C4 (the branches of these solutions are shown in Fig. [1} : C3 for 
Ml = 0.135 (left), C4 for /ii = 0.122 (middle), and C4 for ^ii = 0.123 (right). 
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FIG. 8: (Color online) Top panel: the norm of the numerical (left) and approximate (right) stationary solutions to Eq.([T]) with 
the self-repulsive nonlinearity (cr = 1) as a function of /ii for ^2 = 0.18. Bottom panel: blowups of segments of the top left 
panel where bifurcations of two-components solutions happen. 



14 




FIG. 9; (Color online) The norm of the numerical (left) and semi-analytical (right) solutions of Eq. ([T]) for the self-repulsive 
nonlinearity {a = 1) as a function of ^i, with fi2 = 0.16 (top), /12 ~ 0.14 (middle), ^2 = 0.12 (bottom). The notation is the 
same as in Fig. |8l 
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FIG. 10: (Color online) Profiles of wave functions u and v (solid and dash-dotted lines, respectively, in the top panel), and the 
respective stability eigenvalues (bottom) corresponding to solutions of types CI at ^2 = 0.276 (left), and ASl at /i2 = 0.25 
(right). Branches CI and ASl are shown in Fig. |8l The existence of the quartet of complex stability eigenvalues implies the 
oscillatory instability of the solution. 
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FIG. 11: (Color online) Spatiotemporal contour plots of the densities, \uf and \v\'^, of unstable two-components (combined) 
solutions for the self-attractive nonlinearity (cr = — 1). The left and right panels depict the simulated evolution of wave functions 
It (top) and V (bottom) in the unstable solutions of types C2 and C4, respectively, from Fig. [1] 




FIG. 12; (Color online) The same as in Fig. 1121 but for the self-repulsive nonlinearity (cr = 1). Left and right panels show the 
evolution in the case of the unstable solutions of types C2 and C4, respectively, from Fig. |8] 



